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Abstract. We are interested in a class of numerical schemes for the 
optimization of nonlinear hyperbolic partial differential equations. We 
present continuous and discretized relaxation schemes for scalar, one- 
conservation laws. We present numerical results on tracking type prob- 
lems with nonsmooth desired states and convergence results for higher- 
order spatial and temporal discretization schemes. 
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1 Introduction 

We consider an optimal control problem for scalar conservation laws of the type 

minimize Uo J(u(T),uo) 
subject to u t + f(u) x — 0, u(0,x) = Uq(x), 

Here, J and / are assumed to be smooth and possibly nonlinear functions. The 
initial value uq acts as control to the problem. It can be observed that the wave 
interactions that occur in the solution u in the case of a nonlinear flux function 
/ pose the serious analytical challenges. Recently, the differentiability of J with 
respect to uq could be proven in the sense of shift-differentiability. We refer to 
[6I9I10I11I28I29I30I4I32] for more details. 

Here, a class of numerical methods applied to the optimal control problem 
^ is studied. We only consider the case of smooth initial data and smooth 
solutions u and refer to [4] for more details. For a numerical analysis including 
shock waves and in the case of the Lax-Friedrichs scheme we refer to [2111321 and 
the references therein. 
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1.1 Relaxation Method 

As motiviation for a numerical scheme we follow the ideas of Jin and Xin [22 . 
Therein, a linear approximation (J2J of the nonlinear hyperbolic equation 

d t u + d x f(u) = 

has been discussed. For initial conditions u(x, 0) = uq the approximation is 

d t u + d x v = 0, u(x,0) = uq, 

d t v + a 2 d x u = \ ( f(u) - v) , v(x, 0) = f(u ) (2) 

where e > is the relaxation rate and a is a given constant satisfying the 
subcharacteristic condition max„ < a. For e being small, the solution u 

of © satisfies d t u + d x f(u) = ed x ((a 2 - f{u) 2 )d x u) (cf. [22]). Applying the 
relaxation to the optimal control problem (JT|), we obtain 

{u t + v x = 0, 
Vt + a 2 u x — i (f(u) — v) , (3) 
u(0, x) = u , v(0, x) = f(u ) 

The corresponding adjoint equations for (j3]) are given by (cf. [?]) 

-Vt - a 2 q x = -f'(u), p(T, x) = p T {x), 
e 

-qt-Px = ~, q{T,x) = q T {x). 

e 

For more information on the relaxation system, its limiting scheme for e = 
0, further numerical analysis and extensions we refer to |1I2I5I14 3 8 22 25 27 
and the references therein. Also, the computations are valid provided that all 
appearing functions are at least once differentiable. This is in general not the 
case for conservation laws. 



2 IMEX-Runge-Kutta Discretization 

Numerical discretization of the relaxation system using higher order temporal 
discretizations combined with higher order spatial discretization has been in- 
vestigated in several recent publications as for example [22127] . We apply so 
called implicit-explicit Runge-Kutta methods 26 27J3] as temporal discretiza- 
tion (IMEX RK). Here, the expliciti integration is used for the linear hyperbolic 
transport part and an implicit method is applied to the the stiff source term. 
Implicit-explicit Runge-Kutta method have been studied in the context of con- 
trol problems for example in |4I19[ . Define 



y=(u,v) T , g(y) = (v,a 2 u) T and r(y) := (0, -(v - f(u))f 
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then @ becomes 

Vt + <?(y)x = -r(y), and y(0, x) = (u°, /(^ )) T (a:) 
e 

Applying a suitable discretization D x of the spatial derivative yields the semi- 
discrete state equations 

y' = -D x g(y) + ~r(y), y(0) = y°. (4) 

Remark 1. Spatial discretizations for the linear transport part are well-known. 
The simplest possible is a first-order Upwind method: 

h j = ~i ( ° 2 o) te+va - yi-va) + ;Ky.) , 

where y^+i/2 is obtained by applying the first-order upwind method to char- 
acteristic variables v ± au. Higher order MUSCL schemes, WENO schemes or 
central schemes have also been studied in this context. 

The resulting semi-discrete optimal control problem is then given by: 
minimize j (y(T) , y° ) 

(5) 

subject to y' - -D x g(y) + ±r(y), y(0) = y°. t e [0, T] 

In the context of relaxation schemes the semi-discrete problem is seen as a 
time-integration problem with stiff source which is discretized by an IMEX RK 
methods. For the numerical discretization we therefore consider the previous 
problem as an optimal control problem involving ordinary differential equa- 
tions. Literature concerning the numerical analysis of Runge-Kutta methods 
for the optimality system of ([5]) have been studied in [1717124] . In [7117] parti- 
tioned Runge-Kutta methods for the optimality system are obtained using the 
discretize-then-optimize approach. The derived partitioned Runge-Kutta meth- 
ods have been analysed with regard to symplecticity and order of convergence. In 
[19 , Herty and Schleper, moreover, analysed the associated adjoint imex Runge- 
Kutta method that one obtains if an explicit method is applied to D x g(y) and a 
(diagonally) implicit method to -r(y). In the following, we will analyse general 
partitioned Runge-Kutta methods using IMEX RK methds. More details can be 
found in {20] . Therein, the following IMEX Runge-Kutta discretization of (H|) is 
studied. 

- y n +hj?r\ ay£> x5 (r n W ) +/i£i=i a l3 \r{Y [ n ]) ) i = 1, ..,s 

(6) 

y n +i = Vn + h YX=i &iD x gOti %) ) + /i£i=i ^ir(r„ W ), n = 0, 1,2, . 
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A nonlinear variable transformation and two intermediate states Kn and K$ 
give the equivalent system 

K$ = D x g (y n + h £* =1 aykP + h £' =1 a^K^ i = l,..,s 

K { : ] = \r (y n + h £' =1 OijkW + h £' =1 dijK^) % = 1, .., a (?) 

Vn+l = Vn + hJ2i=i ^iKn + hj^^i ^iKn , Tl = 0, 1,2, . 

The associated optimality systems for the two previous optimization problems 
then coincide and we refer to [3D] for mor details. It is proven that the adjoint 
schemes are equivalent to 



3=1 3=1 6 



S s 1 

PW =p„ -hY^fcgf^fDvPV) -hJ2Pij -r'(Y^) T P (j) * = 1,-,* (8) 



e 

3=1 3=1 



= P« - /l^^ 5 '(^ W ) T ^P W - fc£« i ir'(l r ») T P« ,1 = 0, 1, ..,N-1 



Here, the coefficients of the Runge-Kutta method ay , ay , /3y and /3y are given 
by 

- Wj- _ Wj _ g _ _ Wj Wj 
ay :— — ~Q>ji, ay :— u)j — — djj, py :— cjj a^, py :— Wj ctji. 



2.1 Properties of Discrete IMEX-RK Optimality System 

For the resulting scheme (JSJ,© order conditions can be stated [5UJ. To this end 
we add a suitable equation for p to the previous system. 

s s 

Pn+i = Pn - fy (Y^ ) T P (l) - h £ Ui 9y (y W ) T P^ . (9) 

i=l i=l 

The full method therefore is a standard additive Runge-Kutta scheme for 

y' = -D x g(y) + -r(y) 
e 

p' = </(y) T A c p + V(yfp 
p' = </(y) T A c p + V(yf P 



If we define 



3=1 



Numerical Methods for the Optimal Control of Scalar Conservation Laws 



5 



7* : = J2j=l a i]i and li '■= Y 

J = l 

s 

3=1 

then holds true. 

Theorem 1. Consider the Runge-Kutta scheme ®), (0), (G?P. This scheme is of 

— First- Order : if (SRK1) is of first order 

— Second- Order : if ( SRK1 ) is of second order 

— Third- Order : if (SRK1) is of third order and either 

s ^ Jt y_ . s 1 

Y ^ =3' Y Ui = 3 ' Y Ui = 3' 

i— 1 i— 1 i—1 

are satisfied or if 

s 1 s 1 

Oij 7i = - , 22 U i S *J 7i = g 

t=l i=l 

and i/ 

5 1 s 1 

2 Ui tHj = j or Y u i ~ a v H = c 



6 

<=i 



are satisfied. 



Note that the system ([6]) and (O is not completely coupled, since the forward 
scheme © is solved independently of the adjoint scheme ((5]). General order 
conditions can be found e.g. in [23]. The proof of Theorem [T] and together with 
more details are discussed in [2TJ1. 



3 Numerical Results 
3.1 Scalar Example 

As a simple example, we use a tracking type functional J(u) together with 
Burgers' equation 

and the desired state Ud at final time T = 2.0, that belongs to the initial con- 
dition Ud(0, x) = | + sin (a;) and we start the optimization with the initial data 
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u s (0, x) = 0.5. Moreover, the spatial interval is given by x £ [0,2tt], As dis- 
cretization of the objective functional, we use 



J(u(-,T),u ,u d ) 



Ax 



K 



Udi\ 



Moreover, the discrete gradient of the reduced cost functional is given by 

^u ,J = Po,i + {Df(u ) T qo)i- 

In order to solve the optimal control problem, we apply a steepest descent 
method (with respect to the reduced cost functional) with fixed stepsize < 
a < 1, i.e. we set u^ +1 = Mq + aV„ i J. As stopping criterion for the optimiza- 
tion process we test \J(uo,Ud)\ < tol where tol = IE — 2 denotes a predefined 
stopping tolerance. We observe grid independence in the case where u and u 
are diffcrcntiable in space and time. 

As first-order scheme, we test the Implicit-Explicit Euler scheme 



< +1 =u* -AtD x v* 
< +1 =v* -Ata 2 D x v* 



for the forward, as well as for the backward 



q* = q? +1 - AtD*p? +1 
p * =p n+l_ Ata 2 D * q n 

Q? = it - 4*9? 

vl=vX + f «?/'(«?) 



The spatial gridsize is chosen to be N x = 300, whereas the time discretization is 
done according to the CFL condition with constant ccfl — 0.5. 
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N 


Nr. of It. 


CPU time (in sec.) 


100 


44 


1.795572c+01 


150 


43 


3.768984e+01 


200 


42 


6.380441e+01 


300 


41 


1.491838e+02 



4 Summary 

We briefly discussed a class of numerical methods applied to an optimal control 
problem for scalar, hyperbolic partial differential equations. Order conditions for 
the temporal numerical discretization in the case of differentiable functions have 
been stated. Future work includes the analysis of additional properties of the 
derived numerical discretizations as for example strong stability and asymptotic 
preservation properties. 
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Appendix. 

The discrete adjoint equations that correspond to the discrete optimization prob- 
lem associated with J7]) are 

s s 

3=1 3=1 6 

s s 1 

3=1 3=1 6 

Pn = Pn +1 + £ ffW>) T A$> + E -/(Y^f ■ & 
i=l i=l 

pn =f(yN,y ) ■ 

Moreover, the variable transformation that is needed to obtain © is given by 
and P« := |=- (i=l,..,s; n = 0, .., N - 1) . 
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On the other hand, using ([5]) the associated discrete adjoint equations are 
C« = h (u,if y (YW) + uj t g y (Y^)) T Pn+1 + £ a jt f y (Y^) T C (j) 

3=1 

s 

+hJ2^g v (Y^) T C U) i = l,.., a 

3=1 

s 

Pn=p n+ l+^2Cn\ i=l,..,N-l, PN=j'{VN) 

1=1 

which can be transformed into the scheme ((SJ) using the variable transformation 
PW := + ± ^ and P« := Pn+1 + £ ^ ■ 
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